Sims
Yang Xiao
2024-10-28
Setting
path_proj = here::here()
path_source = file.path(path_proj, "source")
source(file.path(path_source, "simulation", "simulations_functions_final.R"))##
## Attaching package: 'tidyr'## The following object is masked from 'package:reshape2':
##
## smithssource(file.path(path_source, "functions", "plot_function.R"))##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionsource(file.path(path_source, "functions", "fit_function.R"))
source(file.path(path_source, "functions", "table_function.R"))
# place for draws
# mac
posterior_draws_path = file.path(Sys.getenv("HOME"), "Desktop", "draws", "testEach")
#Windows
#posterior_draws_path = file.path(Sys.getenv("USERPROFILE"), "Desktop", "draws", "testEach")
#data path
data_save_path = file.path(path_proj, "data", "fitted_model", "simulation", "4. b_ou")#models
q_constant <- file.path(path_proj, "source", "models",
"q-constant.stan")
b_constant <- file.path(path_proj, "source", "models",
"b-constant.stan")
b_rw <- file.path(path_proj, "source", "models",
"b-rw1.stan")
b_ou <- file.path(path_proj, "source", "models",
"b-ou.stan")
compiled_models <- list(
q_constant = cmdstan_model(q_constant),
b_constant = cmdstan_model(b_constant),
b_rw = cmdstan_model(b_rw),
b_ou = cmdstan_model(b_ou)
)
models_to_use <- c("q_constant", "b_constant", "b_rw", "b_ou")###### setting #####
# data
alpha_increase_seq_1 <- seq(10, 750, by = 30)
alpha_decrease_seq_1 <- seq(750, 10, by = -30)
alpha_lamb = c( rep(10,5), alpha_increase_seq_1 + rnorm(alpha_increase_seq_1,10,10),
alpha_decrease_seq_1 + rnorm(alpha_decrease_seq_1,10,10),
rep(10,5))
beta_lamb = 0.5
T = 60
# reprot delay
D <- 15;
# Time period for checking
D_check <- 5
first_date <- as.Date("2024-01-01")
scoreRange <- c(first_date+days(9), first_date+days(19), first_date+days(29),
first_date+days(39), first_date+days(49))Fully reported case
Simulation
params_b_ou_FR <- list(
data = list(
alpha_lamb = alpha_lamb,
beta_lamb = beta_lamb,
T = T,
date_start = as.Date("2024-01-01"),
D = D
),
q_model = list(
method = "b_ou",
method_params = list(theta_logb = 0.3, mu_logb = log(0.7), init_logb = log(0.7), sigma_logb = 0.2,
theta_logitphi = 0.3, mu_logitphi = 1, init_logitphi = 1, sigma_logitphi = 0.2)
)
)
b_ou_FR <- simulateData(params_b_ou_FR)## [1] "here we go"par(mfrow = c(2, 1))
plot(b_ou_FR$b, pch = 19, type = "b")
plot(b_ou_FR$phi, pch = 19, type = "b")
par(mfrow = c(1, 1))
matplot(t(b_ou_FR$q), type = "l", lty = 1, ylim = c(0, 1))
### Exploratory analysis
# exploritary analysis
page_num <- ceiling(nrow(b_ou_FR$case_reported_cumulated)/16)
exp_plot_ou <- fit_exp_plot(b_ou_FR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# print(exp_plot_ou)
#
# exp_plot_ou$coefficients
exp_b_data_ou<- data.frame( date = as.Date(rownames(b_ou_FR$case_reported_cumulated)),
b = exp_plot_ou$coefficients$b)
exp_b_plot_ou <- ggplot(exp_b_data_ou, aes(x = date, y = b)) +
geom_point(color = "black", size = 1.5) +
geom_smooth(method = "loess", se = TRUE,
color = "blue", fill = "grey", alpha = 0.5) +
theme_minimal() +
labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")
print(exp_b_plot_ou)## `geom_smooth()` using formula = 'y ~ x'
Model fitting
# out_b_ou_FR <-
# nowcasting_moving_window(b_ou_FR$case_reported_cumulated, scoreRange = scoreRange,
# case_true = b_ou_FR$case_true,
# start_date = as.Date("2024-01-01"),
# D = D,
# methods =models_to_use,
# compiled_models = compiled_models,
# iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
# num_chains = 3, suppress_output = T,
# posterior_draws_path = file.path(posterior_draws_path, "b_ou")
# )
#
#
# save(out_b_ou_FR, file = file.path(data_save_path, "FR_b_ou.RData"))
load(file.path(data_save_path, "FR_b_ou.RData"))
results_b_ou_FR <- nowcasts_table(out_b_ou_FR, D = D, report_unit = "day",
methods = models_to_use)
plots_b_ou_FR <- nowcasts_plot(results_b_ou_FR, D = D, report_unit = "day", methods = models_to_use)## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.for (i in 1:length(results_b_ou_FR)) {
print(calculate_metrics(results_b_ou_FR[[i]],
methods = models_to_use))
}## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 85.22 70.55 62.37 69.49 129.22 0.5 q_constant
## 2 19.56 7.57 9.24 5.64 51.91 1.0 b_constant
## 3 23.04 8.33 10.69 5.84 58.30 1.0 b_rw
## 4 19.91 7.60 9.32 5.45 57.60 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 55.60 11.34 38.79 10.77 99.56 0.80 q_constant
## 2 61.74 7.69 25.93 4.09 77.86 0.85 b_constant
## 3 50.65 6.28 21.14 3.30 93.11 0.90 b_rw
## 4 59.77 7.35 23.78 3.63 91.81 0.95 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 207.01 14.54 71.48 6.73 112.68 0.83 q_constant
## 2 190.78 13.20 59.94 4.99 107.87 0.87 b_constant
## 3 24.01 1.84 10.57 1.18 141.98 1.00 b_rw
## 4 54.73 3.82 16.97 1.55 147.04 1.00 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 47.16 4.30 25.56 2.81 114.28 0.95 q_constant
## 2 39.89 3.60 19.18 2.06 111.69 0.95 b_constant
## 3 25.96 2.42 9.50 1.08 132.64 1.00 b_rw
## 4 26.00 2.38 8.43 1.00 131.17 1.00 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 17.38 2.37 12.72 1.83 107.65 1 q_constant
## 2 13.68 2.16 9.82 1.48 106.45 1 b_constant
## 3 12.17 3.01 4.85 1.06 111.39 1 b_rw
## 4 9.97 2.48 4.14 0.94 110.63 1 b_oufor (i in 1:length(results_b_ou_FR)) {
print(calculate_metrics(data.table::last(results_b_ou_FR[[i]],D_check),
methods = models_to_use))
}## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 119.21 69.53 108.52 68.11 222.22 0.2 q_constant
## 2 27.65 10.03 17.66 7.59 84.82 1.0 b_constant
## 3 32.57 11.38 20.69 8.63 96.81 1.0 b_rw
## 4 28.14 10.18 17.88 7.50 95.60 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 91.84 12.78 79.53 11.12 195.01 0.6 q_constant
## 2 123.18 15.02 94.31 11.86 149.82 0.4 b_constant
## 3 101.17 12.37 78.88 9.89 207.62 0.6 b_rw
## 4 119.45 14.53 89.65 11.15 202.43 0.8 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 504.83 34.74 364.40 25.57 262.63 0.2 q_constant
## 2 466.25 32.02 319.70 22.28 246.81 0.2 b_constant
## 3 57.57 4.04 49.10 3.53 444.63 1.0 b_rw
## 4 133.58 9.17 89.63 6.27 472.82 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 115.26 10.27 90.21 8.12 199.41 0.6 q_constant
## 2 98.66 8.76 67.85 6.10 193.81 0.8 b_constant
## 3 71.91 6.56 52.79 4.85 350.24 1.0 b_rw
## 4 72.29 6.44 48.32 4.36 337.03 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 22.96 4.86 20.81 4.26 122.62 1 q_constant
## 2 22.97 5.39 20.94 4.52 119.81 1 b_constant
## 3 37.28 9.35 27.93 6.63 166.41 1 b_rw
## 4 30.57 7.61 24.48 5.68 158.41 1 b_ouprint(plots_b_ou_FR)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
### Parameter Checking
# try the third one, "2024-01-30"
T_test = T * 3/6
# b_ou
varnames_b_ou <- out_b_ou_FR$b_ou[[3]]$summary()$variable
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("sigma_log_b"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("sigma_logit_phi"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("theta_log_b"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("theta_logit_phi"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("mu_log_b"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("mu_logit_phi"), prob = 0.95, prob_outer = 0.95)
param_true = tibble(
parameter = grep("^b\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$b[1:T_test]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("b"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^phi\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$phi[1:T_test]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("phi"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$lambda[1:T_test]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$q[10,]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE)), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_ou_FR$b_ou[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
Not fully reported case
Simulation
# ou_NFR
params_b_ou_NFR <- list(
data = list(
alpha_lamb = alpha_lamb,
beta_lamb = beta_lamb,
T = T,
date_start = as.Date("2024-01-01"),
D = D
),
q_model = list(
method = "b_ou",
method_params = list(theta_logb = 0.2, mu_logb = log(0.1), init_logb = log(0.1), sigma_logb = 0.15,
theta_logitphi = 0.2, mu_logitphi = 1.5, init_logitphi = 1.5, sigma_logitphi = 0.15)
)
)
b_ou_NFR <- simulateData(params_b_ou_NFR)## [1] "here we go"par(mfrow = c(2, 1))
plot(b_ou_NFR$b, pch = 19, type = "b")
plot(b_ou_NFR$phi, pch = 19, type = "b")
par(mfrow = c(1, 1))
matplot(t(b_ou_NFR$q), type = "l", lty = 1, ylim = c(0, 1))
Exploratory analysis
# exploritary analysis
page_num <- ceiling(nrow(b_ou_NFR$case_reported_cumulated)/16)
exp_plot_ou <- fit_exp_plot(b_ou_NFR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
# print(exp_plot_ou)
#
# exp_plot_ou$coefficients
exp_b_data_ou<- data.frame( date = as.Date(rownames(b_ou_NFR$case_reported_cumulated)),
b = exp_plot_ou$coefficients$b)
exp_b_plot_ou <- ggplot(exp_b_data_ou, aes(x = date, y = b)) +
geom_point(color = "black", size = 1.5) +
geom_smooth(method = "loess", se = TRUE,
color = "blue", fill = "grey", alpha = 0.5) +
theme_minimal() +
labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")
print(exp_b_plot_ou)## `geom_smooth()` using formula = 'y ~ x'
Model fitting
# out_b_ou_NFR <-
# nowcasting_moving_window(b_ou_NFR$case_reported_cumulated, scoreRange = scoreRange,
# case_true = b_ou_NFR$case_true,
# start_date = as.Date("2024-01-01"),
# D = D,
# methods =models_to_use,
# compiled_models = compiled_models,
# iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
# num_chains = 3, suppress_output = T,
# posterior_draws_path = file.path(posterior_draws_path, "b_ou")
# )
#
#
# save(out_b_ou_NFR, file = file.path(data_save_path, "NFR_b_ou.RData"))
load(file.path(data_save_path, "NFR_b_ou.RData"))
results_b_ou_NFR <- nowcasts_table(out_b_ou_NFR, D = D, report_unit = "day",
methods = models_to_use)
plots_b_ou_NFR <- nowcasts_plot(results_b_ou_NFR, D = D, report_unit = "day", methods = models_to_use)
for (i in 1:length(results_b_ou_NFR)) {
print(calculate_metrics(results_b_ou_NFR[[i]],
methods = models_to_use))
}## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 12.54 27.74 10.36 22.44 83.51 1.0 q_constant
## 2 63.09 50.82 41.74 48.92 44.51 0.3 b_constant
## 3 65.19 50.87 42.71 49.00 47.51 0.3 b_rw
## 4 65.18 51.47 43.00 49.72 46.60 0.3 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 97.53 20.11 59.79 16.87 105.21 0.75 q_constant
## 2 107.32 21.24 68.16 18.51 98.11 0.70 b_constant
## 3 101.82 19.10 61.49 15.82 151.86 0.85 b_rw
## 4 97.64 19.29 59.39 16.16 137.12 0.85 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 230.91 25.77 164.53 23.15 97.91 0.33 q_constant
## 2 151.09 17.19 91.79 13.17 112.18 0.60 b_constant
## 3 125.71 15.59 68.59 10.97 230.08 0.87 b_rw
## 4 149.25 16.58 84.42 11.93 181.58 0.83 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 208.06 22.57 157.32 19.46 105.31 0.32 q_constant
## 2 122.22 14.75 80.34 11.51 124.76 0.65 b_constant
## 3 120.17 14.57 64.99 9.33 246.94 0.95 b_rw
## 4 110.66 13.29 58.07 8.31 246.22 0.95 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 188.98 22.04 152.18 19.94 99.57 0.24 q_constant
## 2 93.19 12.37 66.32 10.09 114.71 0.60 b_constant
## 3 43.68 10.30 31.00 6.52 224.60 1.00 b_rw
## 4 42.72 9.40 30.72 6.10 211.42 1.00 b_oufor (i in 1:length(results_b_ou_NFR)) {
print(calculate_metrics(data.table::last(results_b_ou_NFR[[i]],D_check),
methods = models_to_use))
}## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 16.91 23.72 15.79 17.89 147.01 1.0 q_constant
## 2 88.46 41.99 72.15 40.41 76.00 0.4 b_constant
## 3 91.46 42.97 74.21 41.21 81.82 0.4 b_rw
## 4 91.43 43.46 74.62 41.89 80.00 0.4 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 187.03 23.63 178.25 22.74 210.24 0.2 q_constant
## 2 203.78 25.76 195.49 24.98 195.03 0.2 b_constant
## 3 195.01 24.79 186.64 23.92 350.04 0.4 b_rw
## 4 187.88 23.68 179.15 22.82 303.84 0.4 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 346.27 24.98 334.78 24.28 214.63 0.0 q_constant
## 2 188.77 13.66 164.09 12.00 248.82 0.4 b_constant
## 3 159.53 11.67 129.36 9.58 735.84 1.0 b_rw
## 4 207.92 15.28 174.11 12.86 517.85 0.8 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 356.28 34.65 342.81 32.82 168.42 0.0 q_constant
## 2 237.49 23.66 202.05 19.77 201.00 0.2 b_constant
## 3 308.36 30.28 286.99 27.68 534.86 0.6 b_rw
## 4 288.14 28.41 262.57 25.42 469.82 0.6 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 125.73 25.42 124.44 25.22 116.81 0.0 q_constant
## 2 49.77 10.41 46.79 9.69 139.01 0.8 b_constant
## 3 79.11 16.99 75.93 15.97 386.05 1.0 b_rw
## 4 82.07 17.20 79.81 16.53 276.44 1.0 b_ouprint(plots_b_ou_NFR)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Parameter Checking
# try the third one, "2024-01-30"
T_test = T * 3/6
# b_ou
varnames_b_ou <- out_b_ou_NFR$b_ou[[3]]$summary()$variable
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("sigma_log_b"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("sigma_logit_phi"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("theta_log_b"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("theta_logit_phi"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("mu_log_b"), prob = 0.95, prob_outer = 0.95)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("mu_logit_phi"), prob = 0.95, prob_outer = 0.95)
param_true = tibble(
parameter = grep("^b\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_NFR$b[1:T_test]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("b"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^phi\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$phi[1:T_test]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("phi"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$lambda[1:T_test]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$q[10,]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_ou, value = TRUE)), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_b_ou, value = TRUE),
x = b_ou_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_ou_NFR$b_ou[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)